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First order phase transitions of spin systems* 

Alain BILLOIRE^ 

'^Service de Physique Theorique de Saclay, 
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1^-^ , I review some numerical ways to determine the parameters of systems close to a first order phase transition 

' point: energy and specific heat of the coexisting phases and interface tension. Numerical examples are given for 
I the 2-d q states Potts model. 
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1. INTRODUCTION 

First order phase transitions are easily de- 
scribed in terms of the competition between two 
phases (see e.g. Q). Consider for definiteness 
a temperature driven phase transition between a 
disordered phase and a q-fold degenerate ordered 
phase of free energy fd{f3) and fo{P)- In order 
to use numerical methods, one needs to know the 
behavior of statistical averages in a finite volume, 
close to the transition point. This was more or 
less correctly guessed by several authors |^,^ , and 
put on rigorous footing (at least for a class of 
models including the large q Potts model) in Ref 
where the partition function is shown to be- 
have like 

ZiP,L) = e-^''^/''(«+ge-^'^^<'('') (1) 
+ 0(e-^^)e-''^(«^' ; 6 > 

This formula is correct for periodic boundary 
conditions, fd{P) and foiP) are L independent 
functions, and /(/?) = min{fd{l3), fo{/3)}. Al- 
though there is no analytical continuation of the 
free energy of one phase into the other phase, 
both fd{P) and /o(/3) are shown to be several 
times differentiable at /3t, the (infinite volume) 
transition point where fo — fd- As a consequence 
the finite size behavior of any moment of the en- 
ergy is obtained by differentiating Eq. |^. For 
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example, the maximum of the energy fluctuation 
behaves like 



with similar results for < S > 



< S >^ 



are simple functions of E^, Ej^, Co and C, 



BL 
lant) 



1^ 



(the Binder Landau cumu- 

. . The constants A and B in Eqj^ 

the 

(infinite volume limit) energy and specific heat of 
the pure phases at (3t. Identical results for the 
first two terms in the 1 /L'^ expansion can be ob- 
tained using the Binder Landau ansatz||^ for the 
energy probability density Pl [E) . 
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l'^PUe-Eq-CqSt)^ 

2^ (3) 

L^pf{E-Ed-CdST)^ 



2Cd 



which can be obtained from Eqj^ by inverse 
Laplace transform, neglecting the 0{e~''^) terms. 
In order to estimate A and B from numerical 
data, the presence of corrections in Eqj^ 

is a nuisance. As first noticed in Q, averages es- 
timated at Pt are better behaved, for example the 
energy fluctuation 



{<E^>- <E >2), 
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A + BL-'^ 



+ 0(6-"^) 



(4) 
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The energy itself behaves hke, 

<E>L ilSt) E >oo {I3t)+0{e-''^) ; b > 0(5) 

A simple powerful method to determine the 
value for Pt follows perform simulations on 
increasing sized lattices, at a value close to f]f 
Extrapolate to neighboring values of /3 using the 
Swendsen Wang reweighting method]^. Pt is de- 
termined as the fixed point of the set of curves 
{< E >L (/3)}- This is similar to Binder's cu- 
mulant method 0| to determine a second order 
phase transition point using the fourth order cu- 
mulant: 1 — ^"^giyi ■ However in the current 
case, the estimate has only exponentially small 
corrections. Application [g[ of this method to 
the 2-d q = 10 Potts model can be found in 
FigJ^. This is a strong first order phase tran- 
sition (see Tab.|l|), the estimated value for /3t is 
1.42596(7) in excellent agreement with the exact 
value ln(vTO-K 1) ~ 1.42606. 
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Co 


Ed — Eo\ 




Cd — Co 






5 
7 
10 
20 


71.3(10) 
18.06(4) 
5.362(3) 


0.05292 
0.35328 
0.69605 
1.19416 


0.0326 
0.2234 
.44763 
.77139 


2512 
48.09 
10.56 

2.70 



Table 1 

The 2-d Potts model: Values of Co from the large 
q expansion, exact values for Ed ~ Eg, Cd — Co, 
and correlation length ^(/3 — > f3^) of the disor- 
dered phase. The correlation length of the or- 
dered phase is smaller. 



Comparison |^ of the predictions of Eq.^ and 
1^ with numerical data from lattices of size 16 < 
L < 100 can be found in Fig.| and FigJ| for the 2- 
d g = 10 Potts model. In both figures the value of 
the specific heat CV divided by 1/L^ is plotted as 
function of l/L"^. The value at L = oo is exactly 
known and the slope is related to the value of Co 
and Cd- An estimate using the large q expansion 
of the model is[|l| Co ^ 18.06(4). (The difference 
Co — Cd is exactly known). The numerical data 
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Figure 1. Estimate of the infinite volume 
transition temperature of the q — 10 Potts 
model obtained from the equality of the in- 
ternal energy of the following pairs of lattices: 
{12, 16},{16, 20},{20, 24},{24, 36},{36, 44},{44, 50}, 
as function of 1/L^, where L is the smallest value 
of each pair. The curve drawn is a fit to the 
shape a + be~'^^. Note that the (unfounded) 
shape a + b/L'^ would give a good fit of the data 
but an incorrect estimate for Pf 



q Co reference 


7 47.5 (25) 1 
50. (10) i 
44. (22) i 




10 10.6 (11) 1 
12.7 (3) i 
~ 18 i 
18.0 (2) i 




15 8.04 (4) 1 


§ 


20 5.2 ( 2) 1 
5.38 ( 4) i 


12 
13 





Table 2 

Numerical estimates of the specific heat of the 
ordered phase of the q states Potts model. Results 
are given in chronological order for each q value. 
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Figure 2. The maximum of the specific heat di- 
vided by , as a function of for the g = 10 
2-d Potts model. The three curves are drawn us- 
ing the values Co = 18.06 x {.95, 1.00, 1.05}. 
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Figure 3. CV{(it)/L'^, as a function of 1/L^ 
for the q — IQ 2-d Potts model. The three 
curves are drawn using the values Co = 18.06 x 
{.95,1.00,1.05}. 



for CV{fit) would allow to estimate correctly the 
value for Co- Using the data for CVmax or BLmin 
would however led to underestimate Co |^,|| • For 
a very strong transition, like the one of the g = 20 
model, results (from similar sized lattices) are un- 
ambiguous ||l||: both Cy(/3t), CVmax and BL,nin 
give estimates for Co that agree with the large q 
result. On the other hand results for the q = 7 
model on lattices as large as L = 64 do not enter 
the asymptotic regime 

Published results for Co from such finite size 
analysis can be found in Tab.^ for the q = 20, 
10 and 7 Potts models. The large q expansion 
estimates can be found in Tab.^ 

2. DETERMINING THE ORDER OF A 
PHASE TRANSITION 

It is clear from the above that equations like 
are only valid in the large volume limit 
L >> ^. Marginal cases (the so called weak first 
order phase transitions) where f is very large and 
the latent heat is very small are quite hard to 
discriminate from second order phase transitions, 
using numerical methods. 

One signature for first order phase transition is 
the " two peak signal" , namely 

i) Pl{E) has two peaks (of heights and 
pmax^ with a minimum between them of height 

pmin 

ii) This structure behaves the way it should as 
L grows, namely the distance between the two 
peaks goes to a constant and e.g. p^nax j pmm _^ 
oo. 

The existence of two peaks in itself is not a 
signal of a first order phase transition. It is 
well known that, for example, the spin probabil- 
ity density of the Ising model has two peaks at 
the (infinite volume) transition temperature [p^ . 
The distance between the peaks goes to zero as 
L grows as L-P/" (here /3 is a critical exponent) 
and the ratio of the peak height to the minimum 
height stays constant. 

The two peak signal is nearly as old as nu- 
merical simulations, it was made more precise 
in Ref.|l|: Define AF = Zn(P™°^/P™™), mea- 
sured at the transition point (at zero field for field 
driven transitions or at the transition tempera- 
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ture for temperature driven transition.). If AF 
increases as L grows the transition is first order, 
otherwise it is second order. 

In the first order case AF behaves as L''"^ for 
the following reason |0,|l9|: Two phase configura- 
tions with almost planar interfaces dominate in 
the region around the minimum of Pl{E). Let 
E — xEo + (1 — x)Ed, X G [0, 1], the contribution 
of such configurations to the partition function is 



\x) 



-2;L''/3/„(/3) -(l-x)L''/3/d(/3) 



X Ve 



e 

d-l 



(6) 



The first two factors are contributions from 
the pure phases. The factor eyip{—2ao,dL'^~^) 
comes from interface effects, ao.d is the order- 
disorder interface tension, and the exponent p 
can be computed in a continuous approximation 
111). Finally AF - 2ao,dL'^~^ However in cur- 
rent simulations this behavior is only seen for ex- 
tremely strong first order phase transitions, like 
the q = 20 Potts model. To summarize, AF grow- 
ing with 1/ is a signal for first order phase transi- 
tion, AF oc L'''~^ is only seen for obviously first 
order cases. 

It is a widespread belief that the behavior of 
BL^in is a better signal for first order phase tran- 
sitions that the behavior of Pl{E) itself. BL^in 
has a non zero limit for first order transitions, and 
a zero limit for second order phase transitions or 
out of the precise transition point. Energy and 
spin fluctuations have similar properties. I would 
like to warn against this belief: 

i) Suppose that, in the simulated range of lat- 
tice size, Pl{E) is approximately L independent, 
BLmin is then a non zero constant. Would any 
one conclude that the transition is first order? 
The only sensible conclusion is rather that one 
does not known the nature of the transition, and 
needs simulations on larger systems, in order to 
see wether Pl (E) will shrink like for a second or- 
der transition, or will develop a two peak struc- 
ture. Note that there do are examples in the lit- 
erature of the use of BLmin to show that a transi- 
tion is first order in cases where Pl{E) is singled 
peaked. 

ii) An example of the use of BLmin to deter- 



mine the order of the 2d q G [3,5] Potts model 
phase transitions can be found in Ref . ||2l[] on lat- 
tices as large as F = 256. BLmin does seem to go 
to a non zero limiting value for q = 5 and to a zero 
limiting value for q = 3 and 4. The exact limit- 
ing value for g = 5 is BLmin = -446 10"'^ (called 
VL,min - I in Fig.4 of Ref.|2^). Looking at the 
data, it means that there must be a crossover for 
L > 256. The success of this method should thus 
be considered as somewhat accidental. 

3. INTERFACE TENSION 

In P^(E) 




Figure 4. Plot of In PL{(3t, E) for selected lattice 
sizes in the q=20 Potts model. The unfolding of 
a fiat region with increasing lattice size L > 38 is 
observed. 



The best known technique to determine the 
order-disorder interface tension is the histogram 
method first introduced in It uses the two 

peak structure of Pl{E), ao.d is estimated as 



2o-„ 



2F' 



1^ 1 U^/^t 



(7) 



with 0{1/L 



rrcctions. This method be- 



came a winner with the introduction of the mul- 
ticanonical algorithm [E3| . This algorithm brings 
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two vast improvements: i) The exponentially 
large slowing down of the Metropolis algorithm 
is strongly reduced, presumably down to a power 
law slowing down, ii) The regions close to the 
peaks and the region close to the minimum are 
sampled equally. Fig.^ shows our results for 
the q = 20 model energy probability distribution. 
One sees the unfolding of a flat region between 
the two peaks, in accordance with expectations 
(Eqj^ in x independent at /3 = /3t). This unfold- 
ing has not (yet) been seen for lower values of 
q. Results for (Tq obtained with this method are 
compared in Tab.^ with the exact result |^ for 
various values of q. Also shown is an estimate 
using measurements of the auto-correlation time 
of the Swendsen-Wang algorithm |^ , which be- 
haves like r ~ Li g'^'^o.dL''-^ rpj^jg estimate 
is in qualitative agreement with the exact result. 
The following two comments are in order 

i) It is somewhat surprising to obtain excellent 
agreement with the exact results for q = 10 and 
q = 7, since in both cases Pl{E) does not unfolds 
its asymptotic shape. 

ii) Ingenious methods |^,^ where two inter- 
faces are created by setting one half of the sys- 
tem at a temperature (3t -\- 5 (it and the other half 
at a temperature I3t — 5(3t give results in strong 
disagreement with the exact result. 
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15 
20 



2cro,d 



Reference 



.1886 (12) 
.20 

.0241 (10) 
.020792 
.094701 
.0978 (8) 
.0950 (5) 
.239234 
.263 (9) 
.370988 
.3714 (13) 



|28| 

|29| 

|14| 

exact 

exact 

lii 

|24| 
exact 

exact 



Table 3 

Comparison of the numerical estimates of the or- 
der disorder interface tension of the 2-d Potts 
model with exact results. 



q=20 Potts model (^=^^) 



003 ,004 



Figure 5. Locations of the maxima of the disor- 
dered peak E'^"'^''^ at (3t in the q=20 Potts model 
as function of the inverse volume. 



Among the by products of interface tension 
measurements are precise determinations of the 
locations of the peaks of Pl{E). These are 
seen to move with L, an effect not described by 
Eqj^. Two models have been proposed to ex- 
plain this dependence. The first one |]l6| involves 
mixed phase effects, the prediction is E^^^'-^{L) ^ 
E^d^i'X)) + A/L. The second one involves higher 
order l/L"^ terms in the Laplace transform of 
Eqj^, neglecting all 0{e^^^) terms, namely all 
mixed phase effects. The result is E"'-^^{L) 
E^^''{oo)+B/L'^, where B is related to the third 
derivative of the free energy of the phase. Numer- 
ical estimate H of E^^-^'^L) as function of l/L'^ 
are shown in Fig js] for the g = 20 model. The line 
drawn is the parameter free prediction of model 
two from Ref. 10 in leading order. Similar results 
are obtained for the q = 10 case. The predicted 
behavior is barely seen, but i) It has been 
shown in Ref. that the discrepancy can be ac- 
commodated by higher orders terms (still without 
any mixed phase effects) ii) A fit to the 1/L form 
gives an estimates of E'^"'^{oo) incompatible with 
the exact value of Ed- 
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